home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / fv_test.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  154 lines
  1. ;$Id: fv_test.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       FV_TEST
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the F-statistic and the probability that two 
  11. ;       vectors of sampled data have significantly different variances. This
  12. ;       type of test is often refered to as the F-variances Test.
  13. ;
  14. ; CATEGORY:
  15. ;       Statistics.
  16. ;
  17. ; CALLING SEQUENCE:
  18. ;       Result = FV_TEST(X, Y)
  19. ;
  20. ; INPUTS:
  21. ;       X:    An n-element vector of type integer, float or double.
  22. ;
  23. ;       Y:    An m-element vector of type integer, float or double.
  24. ;
  25. ; EXAMPLE
  26. ;       Define two n-element vectors of tabulated data.
  27. ;         X = [257, 208, 296, 324, 240, 246, 267, 311, 324, 323, 263, 305, $
  28. ;               270, 260, 251, 275, 288, 242, 304, 267]
  29. ;         Y = [201, 56, 185, 221, 165, 161, 182, 239, 278, 243, 197, 271, $
  30. ;               214, 216, 175, 192, 208, 150, 281, 196]
  31. ;       Compute the F-statistic (of X and Y) and its significance. 
  32. ;       The result should be the two-element vector [2.48578, 0.0540116], 
  33. ;       indicating that X and Y have significantly different variances.
  34. ;         result = fv_test(X, Y)
  35. ;
  36. ; PROCEDURE:
  37. ;       FV_TEST computes the F-statistic of X and Y as the ratio of variances
  38. ;       and its significance. X and Y may be of different lengths. The result 
  39. ;       is a two-element vector containing the F-statistic and its 
  40. ;       significance. The significance is a value in the interval [0.0, 1.0];
  41. ;       a small value (0.05 or 0.01) indicates that X and Y have significantly
  42. ;       different variances.
  43. ;
  44. ; REFERENCE:
  45. ;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
  46. ;       Cambridge University Press
  47. ;       ISBN 0-521-43108-5
  48. ;
  49. ; MODIFICATION HISTORY:
  50. ;       Written by:  GGS, RSI, Aug 1994
  51. ;                    FV_TEST is based on the routine: ftest.c described in 
  52. ;                    section 14.2 of Numerical Recipes, The Art of Scientific 
  53. ;                    Computing (Second Edition), and is used by permission.
  54. ;-
  55.  
  56. function betacf, a, b, x
  57.   on_error, 2
  58.   eps   = 3.0e-7
  59.   fpmin = 1.0e-30
  60.   maxit = 100
  61.   qab = a + b
  62.   qap = a + 1.0
  63.   qam = a - 1.0
  64.     c = 1.0
  65.     d = 1.0 - qab * x / qap
  66.   if(abs(d) lt fpmin) then d = fpmin
  67.   d = 1.0 / d
  68.   h = d
  69.   for m = 1, maxit do begin
  70.     m2 = 2 * m
  71.     aa = m * (b - m) * x / ((qam + m2) * (a + m2))
  72.      d = 1.0 + aa*d
  73.      if(abs(d) lt fpmin) then d = fpmin
  74.      c = 1.0 + aa / c
  75.      if(abs(c) lt fpmin) then c = fpmin
  76.      d = 1.0 / d
  77.      h = h * d * c
  78.      aa = -(a + m) *(qab + m) * x/((a + m2) * (qap + m2))
  79.      d = 1.0 + aa * d
  80.      if(abs(d) lt fpmin) then d = fpmin
  81.      c = 1.0 + aa / c
  82.      if(abs(c) lt fpmin) then c = fpmin
  83.      d = 1.0 / d
  84.      del = d * c
  85.      h = h * del
  86.      if(abs(del - 1.0) lt eps) then return, h
  87.   endfor
  88.   message, 'Failed to converge within given parameters.'
  89. end
  90.  
  91. function gammln, xx
  92.   coff = [76.18009172947146d0,   -86.50532032941677d0,  $
  93.           24.01409824083091d0,    -1.231739572450155d0, $
  94.            0.1208650973866179d-2, -0.5395239384953d-5]
  95.   stp = 2.5066282746310005d0
  96.   x = xx
  97.   y = x
  98.   tmp = x + 5.5d0
  99.   tmp = (x + 0.5d0) * alog(tmp) - tmp
  100.   ser = 1.000000000190015d0
  101.   for j = 0, n_elements(coff)-1 do begin
  102.     y = y + 1.d0
  103.     ser = ser + coff[j] / y
  104.   endfor
  105.   return, tmp + alog(stp * ser / x)
  106. end
  107.  
  108. function ibeta, a, b, x
  109.   on_error, 2
  110.   if (x lt 0 or x gt 1) then message, $
  111.     'x must be in the interval [0, 1].'
  112.   if (x eq 0  or x eq 1) then bt = 0.0 $
  113.   else $
  114.     bt = exp(gammln(a + b) - gammln(a) - gammln(b) + $
  115.              a * alog(x) + b * alog(1.0 - x))
  116.   if(x lt (a + 1.0)/(a + b + 2.0)) then return, $
  117.     bt * betacf(a, b, x) / a $
  118.   else return, 1.0 - bt * betacf(b, a, 1.0-x) / b
  119. end
  120.  
  121. function fv_test, x0, x1
  122.  
  123.   on_error, 2
  124.  
  125.   nx0 = n_elements(x0)
  126.   nx1 = n_elements(x1)
  127.  
  128.   if nx0 le 1 or nx1 le 1 then $
  129.     message, 'x0 and x1 must be vectors of length greater than one.'
  130.  
  131.   type = size(x0)
  132.  
  133.   mv0 = moment(x0)
  134.   mv1 = moment(x1)
  135.  
  136.   if mv0[1] gt mv1[1] then begin
  137.     f = mv0[1] / mv1[1]
  138.     df0 = nx0 - 1
  139.     df1 = nx1 - 1
  140.   endif else begin
  141.     f = mv1[1] / mv0[1]
  142.     df0 = nx1 - 1
  143.     df1 = nx0 - 1
  144.   endelse
  145.  
  146.   prob = 2.0 * ibeta(0.5*df1, 0.5*df0, df1/(df1+df0*f))
  147.  
  148.   if type[2] ne 5 then prob = float(prob)  
  149.     
  150.   if prob gt 1 then return, [f, 2.0 - prob] $
  151.   else return, [f, prob]
  152.  
  153. end
  154.